Quantitative Big Imaging

• Motivation
• Qualitative Approaches
• Thresholding
  • Other types of images
  • Selecting a good threshold
• Implementation
• Morphology
• Contouring / Mask Creation

Applications

• Simple two-phase materials (bone, cells, etc)
  • Beyond 1 channel of depth
• Multiple phase materials
• Filling holes in materials
• Segmenting Fossils
• Attempting to segment the cortex in brain imaging

• Jean Claude, Morphometry with R
  • Online through ETHZ
  • Buy it
• John C. Russ, “The Image Processing Handbook”,(Boca Raton, CRC Press)
  • Available online within domain ethz.ch (or proxy.ethz.ch / public VPN)

Models / ROC Curves

• Julia Evans - Recalling with Precision
• Stripe's Next Top Model

• Exploratory
  • To visually, qualitatively examine samples and differences between them
  • No prior knowledge or expectations
• To test a hypothesis
  • Quantitative assessment coupled with statistical analysis
  • Does temperature affect bubble size?
  • Is this gene important for cell shape and thus mechanosensation in bone?
  • Does higher canal volume make bones weaker?
  • Does the granule shape affect battery life expectancy?

• What we are looking at
• What we get from the imaging modality

• What comes out of our detector / enhancement process

• Identify objects by eye
  • Count, describe qualitatively: “many little cilia on surface”, “long curly flaggelum”, “elongated nuclear structure”
• Morphometrics
  • Trace the outline of the object (or sub-structures)
  • Can calculate the area by using equal-weight-paper
  • Employing the “cut-and-weigh” method

They match up well to the world view / perspective

Model-Based

• \( \rightarrow \) Experimentalist
• Problem-driven
  • Top-down
  • Reality Model-based

Machine Learning Approach

• \( \rightarrow \) Computer Vision / Deep Learning
• Results-driven

• \( F_{system}(a,b) = a*b \)
• \( I_{stimulus} = \textrm{Beam}_{profile} \)
• \( S_{system} = \alpha(\vec{x}) \)

\( \longrightarrow \alpha(\vec{x})=\frac{I_{measured}(\vec{x})}{\textrm{Beam}_{profile}(\vec{x})} \)

In many setups there is un-even illumination caused by incorrectly adjusted equipment and fluctations in power and setups

Frequently there is a fall-off of the beam away from the center (as is the case of a Gaussian beam which frequently shows up for laser systems). This can make extracting detail away from the center challenging

• For absorption/attenuation imaging \( \rightarrow \) Beer-Lambert Law \[ I_{detector} = \underbrace{I_{source}}_{I_{stimulus}}\underbrace{\exp(-\alpha d)}_{S_{sample}} \]

  • Different components have a different \( \alpha \) based on the strength of the interaction between the light and the chemical / nuclear structure of the material \[ I_{sample}(x,y) = I_{source}\exp(-\alpha(x,y) d) \] \[ \alpha = f(N,Z,\sigma,\cdots) \]

• For segmentation this model is:

  • there are 2 (or more) distinct components that make up the image
  • these components are distinguishable by their values (or vectors, colors, tensors, …)

Start out with a simple image of a cross with added noise \[ I(x,y) = f(x,y) \]

The intensity can be described with a probability density function \[ P_f(x,y) \]

By examining the image and probability distribution function, we can deduce that the underyling model is a whitish phase that makes up the cross and the darkish background

Applying the threshold is a deceptively simple operation

\[ I(x,y) = \begin{cases} 1, & f(x,y)\geq0.5 \\ 0, & f(x,y)<0.5 \end{cases} \]

• We can peform the same sort of analysis with this image of cells
• This time we can derive the model from the basic physics of the system
  • The field is illuminated by white light of nearly uniform brightness
  • Cells absorb light causing darker regions to appear in the image
  • Lighter regions have no cells
  • Darker regions have cells

While scalar images are easiest, it is possible for any type of image \[ I(x,y) = \vec{f}(x,y) \]

The intensity can be described with two seperate or a single joint probability density function \[ P_{\vec{f}\cdot \vec{i}}(x,y), P_{\vec{f}\cdot \vec{j}}(x,y) \]

A threshold is now more difficult to apply since there are now two distinct variables to deal with. The standard approach can be applied to both \[ I(x,y) = \begin{cases} 1, & \vec{f}_x(x,y) \geq0.25 \text{ and}\\ & \vec{f}_y(x,y) \geq0.25 \\ 0, & \text{otherwise} \end{cases} \]

This can also be shown on the joint probability distribution as

We can tune the angular acceptance by using the fact \[ \vec{x}\cdot\vec{y}=|\vec{x}| |\vec{y}| \cos(\theta_{x\rightarrow y}) \] \[ I(x,y) = \begin{cases} 1, & \cos^{-1}(\vec{f}(x,y)\cdot \vec{i}) \leq \theta^{\circ} \\ 0, & \text{otherwise} \end{cases} \]

Segmentation and all the steps leading up to it are really a specialized type of learning problem.

Returning to the ring image we had before, we start with our knowledge or ground-truth of the ring

Which we want to identify the from the following image by using image processing tools

What does identify mean?

• Classify the pixels in the ring as Foreground
• Classify the pixels outside of the ring as Background

• True Positive values in the ring that are classified as Foreground
• True Negative values outside the ring that are classified as Background
• False Positive values outside the ring that are classified as Foreground
• False Negative values in the ring that are classified as Background

We can then apply a threshold to the image to determine the number of points in each category

Try a number of different threshold values on the image and compare them to the original classification

• Recall (sensitivity)= \( TP/(TP+FN) \)
• Precision = \( TP/(TP+FP) \)

ROC Curve

Reciever Operating Characteristic (first developed for WW2 soldiers detecting objects in battlefields using radar). The ideal is the top-right (identify everything and miss nothing)

We can then use this ROC curve to compare different filters (or even entire workflows), if the area is higher the approach is better.

Different approaches can be compared by area under the curve

Another way of showing the ROC curve (more common for machine learning rather than medical diagnosis) is using the True positive rate and False positive rate

• True Positive Rate (recall)= \( TP/(TP+FN) \)
• False Positive Rate = \( FP/(FP+TN) \)

These show very similar information with the major difference being the goal is to be in the upper left-hand corner. Additionally random guesses can be shown as the slope 1 line. Therefore for a system to be useful it must lie above the random line.

A more complicated example is finding the road in an image with cars and buildings

From these images, an expert labeled the calcifications by hand, so we have ground truth data on where they are:

We can perform the same analysis on an image like this one, again using a simple threshold to evalulate how accurately we identify the stripes. Since they are white, we can selectively keep the white components.

Since we know the stripes are on the road, we can add another criteria (just the lower half), we can improve the ROC curve substantially.

Going beyond vector the possibilities for images are limitless. The following shows a tensor plot with the tensor represented as an ellipse. \[ I(x,y) = \hat{f}(x,y) \]

Once the variable count is above 2, individual density functions and a series of cross plots are easier to interpret than some multidimensional density hypervolume.

• Here we have a shale sample measured with X-ray tomography with three different phases inside (clay, rock, and air).
• The model is that because the chemical composition and density of each phase is different they will absorb different amounts of x-rays and appear as different brightnesses in the image

Ideally we would derive 3 values for the thresholds based on a model for the composition of each phase and how much it absorbs, but that is not always possible or practical.

• While there are 3 phases clearly visible in the image, the histogram is less telling (even after being re-scaled).

For this exercise we choose arbitrarily 3 ranges for the different phases and perform visual inspection

The relation can explicitly be written out as \[ I(x) = \begin{cases} \text{Void}, & 0 \leq x \leq 0.2 \\ \text{Clay}, & 0.2 < x \leq 0.5 \\ \text{Rock}, & 0.5 < x \end{cases} \]

Matlab

The simplist is a single threshold in Matlab:

thresh_img = gray_img > thresh

A more complicated threshold:

thresh_img = (gray_img > thresh_a) & (gray_img > thresh_b)

Python (numpy)

thresh_img = gray_img > thresh

Python

thresh_img = map(lambda gray_val: gray_val>thresh,
                gray_img)

Java

boolean[] thresh_img = new boolean[x_size*y_size*z_size];
for(int x=x_min;x<x_max;x++)
  for(int y=y_min;y<y_max;y++)
    for(int z=z_min;z<z_max;z++) {
      int offset=(z*y_size+y)*x_size+x;
      thresh_img[offset]=gray_img[offset]>thresh;
    }

In C/C++

bool* thresh_img = malloc(x_size*y_size*z_size * sizeof (bool));

for(int x=x_min;x<x_max;x++)
  for(int y=y_min;y<y_max;y++)
    for(int z=z_min;z<z_max;z++) {
      int offset=(z*y_size+y)*x_size+x;
      thresh_img[offset]=gray_img[offset]>thresh;
    }

Partial Volume Effect

• The partial volume effect is the name for the effect of discretization on the image into pixels or voxels.
• Surfaces are complicated, voxels are simple boxes which make poor representations
• Many voxels are only partially filled, but only the voxels on the surface
• Removing the first layer alleviates issue

• x -> val

• Shown as a function of radius

Returning to the original image of a cross We can now utilize information from neighborhood voxels to improve the results. These steps are called morphological operations.

Like filtering the assumption behind morphological operations are

• nearby voxels in real images are related / strongly correlated with one another
• noise and imaging artifacts are less spatially correlated.

Therefore these imaging problems can be alleviated by adjusting the balance between local and neighborhood values.

A neighborhood consists of the pixels or voxels which are of sufficient proximity to a given point. There are a number of possible definitions which largely affect the result when it is invoked.

• A large neighborhood performs operations over larger areas / volumes
  • Computationally intensive
  • Can smooth out features
• A small neighborhood performs operations over small areas / volumes
  • Computationally cheaper
  • Struggles with large noise / filling large holes

The neighborhood is important for a large number of image and other (communication, mapping, networking) processing operations:

• filtering
• morphological operations
• component labeling
• distance maps
• image correlation based tracking methods

Formally these have two effectively equivalent, but different definitions.

• Pixels which share a face (red line) with the pixel
• Pixels which share a vertex (yellow dot) with the pixel

• Pixels which are distance 1 from the pixel
• Pixels which are distance \( \sqrt{2} \) from the pixel

In 3D there is now an additional group to consider because of the extra-dimension

• Voxels which share a face (red line) with the voxel (6-adjacent)
• Voxels which share an edge (yellow dot) with the voxel (18-adjacent)
• Voxels which share a vertex (purple dot) with the voxel (26-adjacent)

• Voxels which are distance 1 from the voxel
• Voxels which are distance \( \sqrt{2} \) from the voxel
• Voxels which are distance \( \sqrt{3} \) from the voxel

Erosion

If any of the voxels in the neighborhood are 0/false than the voxel will be set to 0

• Has the effect of peeling the surface layer off of an object

Dilation

If any of the voxels in the neigbhorhood are 1/true then the voxel will be set to 1

• Has the effect of adding a layer onto an object (dunking an strawberry in chocolate, adding a coat of paint to a car)